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Efforts  to  discriminate  buried  unexploded  ordnance  from  harmless  surrounding  clutter  are  often  hampered  by  the 
uncertainty  in  the  number  of  buried  targets  that  produce  a  given  detected  signal.  We  present  a  technique  that  helps 
determine  that  number  with  no  need  for  data  inversion.  The  procedure  is  based  on  the  joint  diagonalization  of  a  set  of 
multistatic  response  (MSR)  matrices  measured  at  different  time  gates  by  a  time-domain  electromagnetic  induction 
sensor.  In  particular,  we  consider  the  Naval  Research  Laboratory's  Time-Domain  Electromagnetic  Multisensor  Towed 
Array  Detection  System  (TEMTADS),  which  consists  of  a  5x5  square  grid  of  concentric  transmitter/receiver  pairs.  The 
diagonalization  process  itself  generalizes  one  of  the  standard  procedures  for  extracting  the  eigenvalues  of  a  single 
matrix;  in  terms  of  execution  time,  it  is  comparable  to  diagonalizing  the  matrices  one  by  one.  We  present  the  method, 
discuss  and  illustrate  its  mathematical  basis  and  physical  meaning,  and  apply  it  to  several  actual  measurements  carried 
out  with  TEMTADS  at  a  test  stand  and  in  the  field  at  the  former  Camp  Butner  in  North  Carolina.  We  find  that  each 
target  in  a  measurement  is  associated  with  a  set  of  nonzero  time-dependent  MSR  eigenvalues  (usually  three),  which 
enables  estimation  of  the  number  of  targets  interrogated.  These  eigenvalues  have  a  characteristic  shape  as  a  function  of 
time  that  does  not  change  with  the  location  and  orientation  of  the  target  relative  to  the  sensor.  We  justify  analytically 
and  empirically  that  symmetric  targets  have  pairs  of  eigenvalues  with  constant  ratios  between  them 
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ABSTRACT 

Efforts  to  discriminate  buried  unexploded  ordnance  from 
harmless  surrounding  clutter  are  often  hampered  by  the  uncer¬ 
tainty  in  the  number  of  buried  targets  that  produce  a  given 
detected  signal.  We  present  a  technique  that  helps  determine 
that  number  with  no  need  for  data  inversion.  The  procedure 
is  based  on  the  joint  diagonalization  of  a  set  of  multistatic  re¬ 
sponse  (MSR)  matrices  measured  at  different  time  gates  by  a 
time-domain  electromagnetic  induction  sensor.  In  particular, 
we  consider  the  Naval  Research  Laboratory’s  Time-Domain 
Electromagnetic  Multisensor  Towed  Array  Detection  System 
(TEMTADS),  which  consists  of  a  5  X  5  square  grid  of  con¬ 
centric  transmitter/receiver  pairs.  The  diagonalization  process 
itself  generalizes  one  of  the  standard  procedures  for  extracting 


the  eigenvalues  of  a  single  matrix;  in  terms  of  execution  time,  it 
is  comparable  to  diagonalizing  the  matrices  one  by  one.  We  pre¬ 
sent  the  method,  discuss  and  illustrate  its  mathematical  basis 
and  physical  meaning,  and  apply  it  to  several  actual  measure¬ 
ments  carried  out  with  TEMTADS  at  a  test  stand  and  in  the  field 
at  the  former  Camp  Butner  in  North  Carolina.  We  find  that  each 
target  in  a  measurement  is  associated  with  a  set  of  nonzero  time- 
dependent  MSR  eigenvalues  (usually  three),  which  enables  es¬ 
timation  of  the  number  of  targets  interrogated.  These  eigenva¬ 
lues  have  a  characteristic  shape  as  a  function  of  time  that 
does  not  change  with  the  location  and  orientation  of  the  target 
relative  to  the  sensor.  We  justify  analytically  and  empirically 
that  symmetric  targets  have  pairs  of  eigenvalues  with  constant 
ratios  between  them. 


INTRODUCTION 

The  United  States  military  has  identified  as  one  of  its  main  en¬ 
vironmental  challenges  reducing  the  worldwide  profusion  of  unex¬ 
ploded  ordnance  (UXO)  contaminating  former  military  training 
grounds  and  battle-scarred  territory.  UXO  kill  or  maim  more  people 
in  conflict  zones  each  year  than  do  landmines  (Moyes  et  al.,  2002); 
they  also  cause  frequent  disruption  in  nations  at  peace.  Remediation 
of  UXO-tainted  land  is  laborious  and  expensive:  targets  of  interest 
are  usually  surrounded  by  metallic  clutter,  natural  and  artificial,  that 
is  also  detected  by  sensors,  and  in  some  areas,  the  targets  them¬ 
selves  are  also  densely  clustered.  Thus,  it  is  often  not  clear  how 
many  buried  targets  contribute  to  a  given  detected  signal,  or  how 


many  are  dangerous,  and  this  shortage  of  information  complicates 
identification  attempts.  This  paper  aims  to  ease  UXO  remediation 
by  presenting  a  method  to  estimate  the  number  of  targets  and  thus 
assess  scene  complexity  quickly  and  reliably  on  the  basis  of  elec¬ 
tromagnetic  induction  (EMI)  data. 

In  EMI  sensing,  a  time-dependent  primary  magnetic  field  estab¬ 
lished  by  a  sensor  induces  eddy  currents  and/or  magnetic  domain 
realignments  in  nearby  ferrous  and  nonferrous  metallic  objects. 
This  response,  in  turn,  gives  rise  to  a  measurable  secondary  mag¬ 
netic  field.  The  measured  signal,  often  referred  to  as  an  “anomaly” 
detected  at  a  “cell,”  depends  on  the  sizes,  shapes,  and  material 
properties  of  the  targets  under  interrogation  and  also  on  the  parti¬ 
culars  of  the  measurement,  i.e.,  the  relative  distances  and  attitudes 
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between  the  targets  and  the  sensor.  The  main  task  of  the  data  ana¬ 
lysis  is  to  separate  these  “intrinsic”  and  “extrinsic”  factors  to  deter¬ 
mine  what  targets  are  producing  the  anomaly  and  where  and  how 
deep  they  are  within  the  cell.  This  inverse  problem  is  usually  solved 
using  a  combination  of  linear  and  nonlinear  least  squares  to  fit  a 
model  to  the  measurements  (Aster  et  al.,  2005;  Everett,  2012). 

Several  models  used  with  success  to  detect  and  characterize  bur¬ 
ied  metal  targets  build  on  the  premise  that  the  currents  and  domains 
that  generate  electromagnetic  response  are  not  distributed  uniformly 
within  each  target,  but  instead,  tend  to  concentrate  around  certain 
singular  points  or  regions.  The  techniques  synthesize  that  response 
as  though  it  consisted  of  a  set  of  analytic  solutions  of  the  Maxwell 
equations  due  to  point  sources  located  at  those  so-called  “scattered 
field  singularities”  (Bogdanov  et  al.,  1999;  Karkashadze  et  al., 
2009).  The  time-dependent  behavior  of  these  sources  is  used  to 
characterize  and  classify  the  causative  targets.  (In  this  work,  we 
use  “source”  to  refer  to  these  elementary  analytic  solutions,  and 
not,  as  is  often  the  case  in  geophysics,  to  transmitters.) 

The  number  of  scattered  field  singularities  depends  on  the  wave¬ 
length  and  on  the  scatterers’  geometry  and  material  properties.  In 
the  EMI  frequency  regime,  where  the  wavelengths  are  significantly 
larger  than  the  targets’  characteristic  lengths,  most  homogeneous 
compact  targets  have  a  dominant  singularity  at  their  centers.  There¬ 
fore,  the  simplest  and  most  widely  used  approach  replaces  each  tar¬ 
get  with  a  single  point  magnetic  dipole  located  near  its  geometric 
center  (Das  et  al.,  1990;  Barrow  and  Nelson,  2001;  Bell  et  al.,  2001 ; 
Pasion  and  Oldenburg,  2001;  Zhang  et  al.,  2003;  Smith  and  Mor¬ 
rison,  2004;  Tarokh  et  al.,  2004;  Fernandez  et  al.,  2011).  This  point 
dipole  model  tends  to  work  best  for  relatively  small  and  distant  tar¬ 
gets,  and  in  terms  of  predictive  accuracy  and  consistency  has  been 
superseded  by  generalizations  that  distribute  dipole-moment  densi¬ 
ties  over  surfaces  or  volumes  (Shubitidze  et  al.,  2010b,  2010c).  On 
the  other  hand,  all  of  these  models  fare  less  well  when  analyzing 
anomalies  due  to  more  than  one  target:  The  inversion  becomes  more 
computationally  expensive,  has  more  local  minima,  and  may  re¬ 
quire  regularization  (Song  et  al.,  2009,  2011;  Grzegorczyk  et  al., 
2011).  When  the  number  of  targets  is  not  known  a  priori,  as  in 
all  real-world  tests,  the  situation  is  complicated  further  because 
the  models  must  be  run  several  times  (assuming  different  target 
numbers)  and  have  the  tendency  to  overfit  the  data. 

The  approach  we  introduce  here  is  designed  to  make  a  fast  es¬ 
timate  of  the  number  of  targets  in  a  cell.  It  is  not  meant  to  replace 
actual  inversion  because  it  does  not  yield  location  information  and 
does  not  characterize  targets  with  complete  precision,  but  rather 
constitutes  an  initial  pre-inversion  step.  (The  inversion  routines  that 
use  as  input  the  resulting  number  estimates  are  presented  elsewhere 
(Shubitidze  et  al.,  2010a,  2010b)  and  will  be  studied  in  detail 
in  a  companion  paper.) 

The  technique  is  based  on  the  joint  diagonalization  (JD)  of  a  se¬ 
quence  of  square  time-dependent  multistatic  response  (MSR)  ma¬ 
trices  that  are  synthesized  directly  from  measured  EMI  signals 
without  invoking  a  forward  model.  The  number  of  nonzero  time- 
dependent  eigenvalues  of  the  set  of  matrices  is  related  to  the  number 
of  meaningful  elementary  sources  present  in  the  illuminated  cell.  As 
will  be  shown  in  the  “Methods”  section,  a  single  source  generating 
three  eigenvalues  may  suffice  to  describe  a  small  or  deeply  buried 
target.  A  complex  field  structure,  caused  by  a  large,  heterogeneous, 
or  very  shallow  target,  may  require  more  than  one  elementary 
source  (and  thus  more  than  three  time-dependent  eigenvalues)  to 


describe  the  resulting  field  adequately.  Objects  like  thin  wires  or 
rings  that  are  “small”  along  one  or  more  dimensions  may  require 
less  than  three  time-dependent  eigenvalues. 

Joint  diagonalization  has  become  an  important  tool  for  signal 
processing  and  inverse  problems,  including  independent  component 
analysis  (Comon,  1994),  blind  source  separation  (BSS)  (Belouchra- 
ni  et  al.,  1997),  common  principal  component  analysis,  and  kernel- 
based  nonlinear  BSS  (Harmeling  et  al.,  2003).  To  the  best  of  our 
knowledge,  this  is  the  first  use  of  this  technique  for  UXO  discrimi¬ 
nation. 

METHODS 

The  TEMTADS  sensor  array 

EMI  sensors  operate  at  very  low  frequencies  (a  few  Hz  to  a  few 
hundred  kHz)  to  penetrate  conductive  ground,  resulting  in  data  with 
poor  spatial  resolution.  Advanced  EMI  sensors  have  been  designed 
to  overcome  this  limitation:  They  enhance  the  information  content 
of  the  data  by  increasing  the  number  of  sensor/target  attitudes  and 
reduce  positioning  uncertainty  by  employing  rigid  arrays  of  trans¬ 
mitters  and  receivers.  One  such  instrument  is  the  Time-domain 
Electromagnetic  Multisensor  Towed  Array  Detection  System 
(TEMTADS),  developed  by  the  Naval  Research  Laboratory  and 
G&G  Sciences,  Inc.  (Steinhurst  et  al.,  2010),  to  which  we  restrict 
attention. 

Shown  in  Figure  1 ,  TEMTADS  is  a  time-domain  EMI  sensor  that 
consists  of  25  transmitter/receiver  (Tx/Rx)  pairs,  each  composed  of 
a  35-cm  square  transmitter  loop  surrounding  a  25-cm  square  recei¬ 
ver  loop.  The  Tx/Rx  pairs  are  arrayed  as  a  rectangular  5x5  grid 
with  40-cm  neighbor-to-neighbor  separation  for  a  total  array  dimen¬ 
sion  of  2  m.  At  each  location,  the  sensor  activates  the  transmitter 
loops  in  sequence,  and  for  each  transmission  all  receivers  record 
signals,  providing  625  EMI  transients  from  ~100  ps  to  25  ms  over 
Ng  =  123  time  gates. 

Multistatic  response  matrices 

To  construct  the  MSR  matrices,  one  assembles  the  625  readings 
at  each  time  gate  into  a  25  X  25  array  so  that  each  column  stands  for 
one  of  Nt  transmitters  and  each  row  represents  one  of  Nr  receivers: 


'  H  „ 

Hn  ■ 

■  ■  H  Wt  - 

H2\ 

H22  ■ 

■  ■  h2Ni 

"1 

_ 1 

HNr  2  ■ 

^NrNt  _ 

where  the  element  Htj  is  the  field  measured  by  the  ith  receiver  in 
response  to  the  jth  transmitter.  Below,  we  look  at  the  time  devel¬ 
opment  of  the  MSR  sequence;  here,  we  concentrate  on  a  single  time 
gate  and  interpret  the  information  contained  in  each  matrix.  It  is 
shown  in  this  paper  that  the  measured  data  are  resolved  as  a  super¬ 
position  of  “elemental”  spatio-temporal  subsignals,  each  related  to 
the  EMI  signature  of  an  elementary  magnetic  dipole  source.  As¬ 
sume  there  are  N  such  sources.  We  postulate  that  the  field  of  the 
jth  transmitter,  immediately  upon  being  shut  off,  induces  in  the 
/th  source  a  dipole  moment  given  by 

m/V  =  U/A;Uf  B*  (2) 


Downloaded  16  Aug  2012  to  129.170.195.144.  Redistribution  subject  to  SEG  license  or  copyright;  see  Terms  of  Use  at  http://segdl.org/ 


Joint  diagonalization  and  UXO 


WB151 


where  Bpzr  represents  the  primary  magnetic  induction  and  the  Euler 
rotation  matrix  U/  relates  the  transmitter  array  coordinate  axes  to  the 
principal  axes  of  the  source.  The  diagonal  polarizability  matrix  A/, 
intrinsic  to  the  source,  measures  the  strength  of  the  moment  that  the 
primary  field  induces  along  each  of  the  source  axes.  This  polariz¬ 
ability  matrix  contains,  in  principle,  all  the  time-dependent  informa¬ 
tion  in  the  measured  secondary  signal. 

To  compute  the  primary  magnetic  induction  Bp;r,  recall  that  the 
TEMTADS  transmitter  assembly  consists  of  coplanar  square  loops 
forming  a  regular  grid.  The  Biot-Savart  law  prescribes  the  primary 
magnetic  induction  at  location  rz  of  the  /th  source  when  the  /th 
transmitter  antenna  (of  area  <7Xx  )  is  excited  by  a  current  If 


■ppr  f*0  h 


Jjl 


1 
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where  pi0  is  the  permeability  of  free  space,  d\'  is  the  usual  Biot- 
Savart  line  element,  and  we  have  introduced  the  normalized  Green 
function  gp;r. 

According  to  Faraday’s  law,  the  signal  measured  by  a  receiver 
coil  is  the  negative  time  derivative  of  the  magnetic  flux  passing 
through  it.  (Long  enough  after  transmitter  shut-off,  this  flux  is 
purely  secondary.)  The  field  at  r  of  a  dipole  of  moment  m  at  rz  is 


u  o  _  /  r  —  r/ 

B  =  — Vx  m  x  - - 
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and  thus 


r-r 


(4) 


by  straightforward  application  of  Stokes’s  theorem.  Thus,  the  signal 
sampled  at  time  tk  by  the  /th  receiver  (of  area  erRx )  when  the  /th 
source  is  excited  by  the  jth  transmitter  is 


^//(^Wx^Tx/y 


Mo  J_  f  4Tx(r'-r,) 
4^ffRx'  (tRXj  ,/RXi  |r'  — r;|3 


m ji(h) 


=  8/^Rx,  •  mjfa) 

=  g^RXi  •  [UA,(4)Ur]  •  g^oTx.Ij,  (5) 


where  a  dot  over  a  variable  indicates  its  time  derivative.  We  have 
also  introduced  the  normalized  Green  function  g|9  in  analogy  with 
g p[;  the  two  are  seen  to  have  the  same  form.  In  equations  4  and  5,  the 
line  element  d\'  lies  on  the  x-y  plane,  and,  as  a  consequence,  the 
Green  functions  are  similar  in  structure  to  those  of  the  simple  model 
presented  in  Appendix  A  (in  which  vertical  dipoles  establish  the 
primary  field  and  only  the  z-component  of  the  secondary  field  is 
measured).  Note  that,  in  the  definition  of  the  received  signal,  we 
have  included  the  exciting  current  Ij  and  the  transmitter  and  recei¬ 
ver  areas;  these  quantities  are  known  and  can  be  factored  out. 

At  a  given  time  gate,  and  with  only  the  /th  source  present,  we 
construct  the  MSR  matrix  for  the  complete  transmitter/receiver  ar¬ 
ray  by  tiling  all  the  Nr  xNt  available  samples  of  expression  5: 

S  =  GscU/A/Uf  (Gpr)r  =  (GscU/)A/(GprU/)r,  (6) 


where  the  primary  (or  transmitter)  dyad  Gpr  is  of  size  Nt  X  3,  the 
secondary  (or  receiver)  dyad  Gsc  is  of  size  Nr  X  3,  and  the  response 
matrix  U/AUf  is  3  X  3.  The  matrix  S  has  size  Nr  X  Nt  and  is  square 


if  Nr  =  Nt,  as  with  TEMTADS.  Moreover,  the  condition 
G^c  =  Gpr  =  G /  holds  approximately  even  though  the  receivers 
do  not  coincide  exactly  with  the  transmitters.  This  means  that 
the  matrix  can  be  diagonalized  to  yield  real  eigenvalues  and 
eigenvectors.  (This  diagonalization  is  more  precisely  a  singular  va¬ 
lue  decomposition  (SVD)  for  each  time  gate,  and  we  use  “diago¬ 
nalization”  as  shorthand  for  “SVD  of  a  symmetric  matrix”  in 
what  follows.)  Because  U/A/Uf  is  3  X  3,  the  matrix  S  has  only  three 
nonzero  eigenvalues  for  a  single  source;  this  property  is  the  basis  of 
the  method.  When  there  is  more  than  one  source  present,  the  de¬ 
composition  6  expands  to 
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where  again  G^c  =  Gpr  =  G /.  Formally,  the  intrinsic  middle  matrix 
in  equation  7  has  size  3 N  X  3 N,  where  N  is  the  number  of  sources; 
numerically,  it  is  a  3V  X  3 N  block-diagonal  matrix  padded  with 
zeros  to  reach  a  size  of  Nr  X  Nt  and  thus  with  only  3 N  nonzero 


Figure  1.  Photograph  (a)  and  sketch  (b)  of  the  TEMTADS  sensor 
array. 
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eigenvalues.  The  method  should  be  able  to  resolve  up  to  [Nr/3\ 
responding  dipole  sources,  or  eight  for  TEMTADS. 

This  procedure  would  be  equivalent  to  full  inversion  (though 
without  the  positional  information)  if  one  could  extract  directly 
the  derivative  A  of  the  polarizability  matrix.  This  is  not  possible, 
however,  because  the  Green  dyads  are  not  orthogonal.  The  actual 
quantities  obtained  from  the  diagonalization  can  be  found  by  repla¬ 
cing  the  Green  dyads  by  their  SYDs: 

S  =  GUAUrGr  =  W[ZVrUAUrVZ]Wr 

=  WZAZrWr  =  YAYr.  (8) 

(The  matrix  within  brackets  in  the  second  step  is  real  and  symmetric 
and  thus  has  a  purely  real  eigendecomposition.)  Equation  8  shows 
that  the  eigenvalues  stored  in  matrix  A  are  not  solely  composed  of 
source  responses,  but  also  contain  location  and  orientation  informa¬ 
tion  from  the  Green  dyads.  The  crucial  point  is  that  the  eigenvalues 
contain  all  of  the  time  dependence  of  the  signal,  a  property  that  we 
now  exploit. 

Time  development  of  the  MSR  eigenvalues 

The  second  step  of  the  procedure  consists  of  incorporating  the 
time  dependence  of  the  MSR  matrices  in  the  eigenvalue  analysis 
from  the  previous  section.  As  we  shall  see  below,  this  information 
is  more  useful  if  each  eigenvalue  can  be  tracked  separately  as  the 
signal  decays  and  associated  with  the  same  time-independent  eigen¬ 
vector.  Thus,  it  is  necessary  to  diagonalize  the  Ng  =  123  matrices 
simultaneously  such  that  they  share  the  same  set  of  orthonormal 
time-independent  eigenvectors,  and  these  appear  in  the  same  order 
at  all  times.  Denoting  again  as  S (tk)  the  MSR  matrix  at  the  kth  time 
gate,  a  unitary  matrix  V  is  sought  such  that  the  products 

Djt  =  WS(4)V  (9) 

are  “as  diagonal  as  possible”  in  the  sense  of  definition  B-l  in  Ap¬ 
pendix  B.  From  now  on,  we  concentrate  on  these  time-decaying 
diagonal  elements,  referring  to  them  as  “eigenvalues”  even  though 


Figure  2.  A  test- stand  experiment  featuring  an  increasing  number 
of  targets  under  interrogation  by  TEMTADS.  As  in  Figure  1,  the 
system  is  observed  essentially  head-on  from  the  +y-direction. 


eigenvalues  are  in  rigor  scalars,  not  curves,  and  describing  them  as 
“parallel”  or  “crossing  each  other”  in  a  way  that  the  figures  will 
make  clear. 

One  could  in  principle  diagonalize  the  MSR  matrix  at  each  time 
gate,  and  the  eigenvectors,  which  depend  only  on  geometry,  should 
stay  constant.  The  fact  that  we  cannot  know  a  priori  the  order  in 
which  the  eigenvalues  and  their  corresponding  eigenvectors  result 
from  each  single-gate  diagonalization  would  be  only  a  minor  com¬ 
plication  if  the  data  had  no  noise.  Instead,  we  look  for  an  orthogonal 
matrix  of  eigenvectors  that  diagonalizes  all  the  MSR  matrices  si¬ 
multaneously.  Two  different  procedures,  described  in  Appendix  B, 
were  employed  with  identical  results.  The  most  efficient  one 
(Flury  and  Gautschi,  1986;  Cardoso  and  Souloumiac,  1996;  Belou- 
chrani  et  al.,  1997)  is  a  well-known  generalization  of  the  Jacobi 
method  for  diagonalizing  single  matrices.  We  note  that  further  de¬ 
velopment  of  our  procedure,  generalized  for  nonsymmetric  and 
nonsquare  MSR  matrices,  may  require  the  use  of  a  joint  SYD  algo¬ 
rithm,  with  two  sets  of  shared  vectors  instead  of  one  (Maehara  and 
Murota,  2011). 

RESULTS 

Test-stand  measurements 

The  main  purpose  of  the  JD  technique  is  to  estimate  the  number 
of  targets  producing  a  given  measured  signal.  In  this  section,  we 
first  describe  an  idealized  experiment  that  was  designed  to  act  as 
a  benchmark  for  multitarget  scenarios  while  avoiding  some  of 
the  complications  that  arise  in  the  field.  The  experiment  was  carried 
out  by  personnel  from  Nova  Research,  Inc.  on  15  January  2010.  The 
setup  is  depicted  in  Figure  2.  The  TEMTADS  sensor  array  was 
placed  on  an  elevated  platform.  Underneath,  several  targets  were 
added  in  sequence  to  form  an  increasingly  complex  scene,  with 
the  instrument  measuring  the  response  of  every  configuration. 

The  targets  of  interest  (TOI)  were  two  common  UXO,  a  105 -mm 
M60  howitzer  shell  and  a  60-mm  mortar  round.  Each  lay  on  its  side 
with  nose  pointing  in  the  +x  direction:  the  105-mm  was  at  a  depth 
of  63  cm  below  the  center  of  the  sensor,  while  the  60-mm  was  38  cm 
below  the  instrument  and  displaced  40  cm  in  the  — y-direction  from 
the  sensor  center.  A  sphere  and  a  spheroid  were  successively  placed 
27.5  cm  below  TEMTADS  to  simulate  shallow  clutter  with  the 
potential  to  obscure  the  targets  of  interest.  The  sphere,  placed  at 
10(x  +  y)  cm  from  the  origin,  was  made  of  aluminum  and  was 
10  cm  in  diameter;  the  spheroid,  made  of  steel  with  major  and  minor 
axes  20  and  4  cm,  respectively,  was  at  — (30x  +  40y)  cm  from  the 
center. 

Initially,  the  targets  were  interrogated  individually.  Figure  3a 
shows  the  complete  set  of  25  eigenvalues  obtained  for  the  105-mm 
shell,  Figure  3b  shows  the  same  result  for  the  60-mm  round,  and 
Figure  4  displays  the  decay  curves  for  the  sphere  and  spheroid 
added  in  the  course  of  the  experiment.  In  all  cases,  at  least  three 
time-dependent  eigenvalues  emerge  from  the  background  noise. 
(Eigenvalues  that  could  be  associated  with  high-order  sources  such 
as  quadrupoles  are  indistinguishable  from  the  noise.)  Two  eigenva¬ 
lues  look  parallel  to  one  another  (i.e.,  there  is  a  constant 
ratio  between  the  two);  for  the  sphere  of  Figure  4a  two  coincide 
and  the  third  is  parallel.  This  parallelism  is  indicative  of  cylindrical 
symmetry,  as  we  justify  in  Appendix  A  using  a  simplified  model. 
The  60-mm  round  has  eigenvalues  of  comparable  size  to  those  of 
the  105-mm  shell,  even  though  the  target  is  much  smaller.  This  is 
due  to  the  fact  that  the  eigenvalues  do  not  store  only  intrinsic 
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responses  but  also  positional  information.  It  can  also  be  noted  that 
the  105 -mm  shell  has  at  least  one  more  above-noise  eigenvalue, 
showing  that  it  is  not  as  accurately  represented  by  a  single  source 
as  the  60-mm  round. 

The  fact  that  the  eigenvalues  change  with  location  points  out  the 
necessity  of  diagonalizing  the  sequence  of  MSR  matrices  simulta¬ 
neously.  Consider  the  60-mm  round  of  Figure  3b.  TEMTADS  mea¬ 
surements  were  carried  out  on  the  same  target  at  a  different  set  of 
locations:  Under  the  sensor  at  a  depth  of  35  cm  (3  cm  shallower  than 
in  Figure  3b),  —  10(x  +  y)  cm  from  the  center  and  29.5  cm  below 
the  sensor,  and  at  the  latter  depth  but  under  the  center.  Figure  5 
shows  the  eigenvalue  decay  curves  that  result  from  carrying  out 
joint  diagonalization  on  these  three  data  sets,  whereas  Figure  6 
shows  the  same  decay  curves  computed  time  gate  by  time  gate  using 
the  SVD.  Both  methods  find  three  above-threshold  eigenvalues  and 
correctly  predict  that  there  is  only  one  target,  but  the  eigenvalues’ 
time  developments  have  a  crucial  difference:  The  JD  eigenvalues 
preserve  their  shapes  case  to  case,  even  as  the  relative  amplitudes 
vary,  whereas  those  obtained  with  the  SYD  change  shape  as  the 


object  location  varies.  The  simple  sorting  provided  by  the  SVD  ex¬ 
plicitly  disallows  the  eigenvalue  crossings  that  help  preserve  their 
shapes  and  the  parallelism  that  reveals  possible  symmetry;  a  more 
systematic  procedure  based  on  eigenvector  comparison  (as  opposed 
to  eigenvalue  comparison)  works  well  for  one-target  scenarios,  but 
breaks  down  in  multitarget  cases  where  there  are  many  more  decay 
curves  crossing  each  other.  Again,  the  number  of  sources  would  be 
estimated  correctly,  but  potentially  important  information  given  by 
the  decay  shape  would  be  lost. 

Figure  7a  shows  the  eigenvalue  decay  curves  for  the  configura¬ 
tion  consisting  of  the  two  TOI  at  the  locations  from  Figure  3.  There 
are  seven  above-noise  eigenvalue  curves,  indicating  that  there  are  at 
least  two  targets  in  the  data.  The  two  sets  of  crossing  decay  curves 
from  Figure  3  are  visible  now  as  a  weighted  nonlinear  superposition 
of  the  individual  responses,  including  the  fourth  eigenvalue  asso¬ 
ciated  with  the  larger  shell.  Figure  7b  and  7c  shows  the  MSR  ei¬ 
genvalues  obtained  after  sequentially  adding  the  first  and  second 
clutter  items.  Ten  eigenvalue  curves  are  above  the  noise  in  Figure  7b, 
revealing  the  presence  of  three  targets,  and  12  curves  appear  in 


Time  (ms)  Tjme  (ms) 


Figure  3.  (a)  Eigenvalues  of  the  TEMTADS  multistatic  response 
matrix  for  a  105 -mm  M60  howitzer  shell  (shown  in  the  inset), 
(b)  Eigenvalues  of  the  TEMTADS  multistatic  response  matrix 
for  a  60-mm  M49  mortar  round  (shown  in  the  inset). 


Figure  4.  Measured  (TEMTADS  multistatic  response)  eigen¬ 
values  for  an  aluminum  sphere  of  diameter  10  cm  (a)  and  a 
steel  prolate  spheroid  (b)  of  minor  and  major  axes  4  and  20  cm, 
respectively. 
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Figure  7c,  showing  that  the  number  has  increased.  The  eigenvalue 
signatures  of  the  original  targets  can  still  be  discerned,  but  they  are 
mixed  with  contributions  from  clutter  items.  A  scenario  with  more 
than  six  nonzero  eigenvalues  potentially  indicates  the  presence  of  a 
single  large  shallow  target  or  of  several  smaller  ones  and  will  be 
routinely  dug  out  when  encountered  in  the  field. 


Camp  Butner  blind  test 

As  part  of  its  effort  to  make  UXO  remediation  faster  and  more 
economical,  the  Environmental  Strategic  Technology  Certification 
Program  (ESTCP)  is  in  the  process  of  administering  a  series  of  blind 
classification  tests  of  increasing  realism  to  measure  the  progress 
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Figure  5.  Eigenvalues  for  the  60-mm  round  computed  using  joint 
diagonalization  for  three  different  object  locations. 


Figure  6.  Eigenvalues  for  the  60-mm  round  computed  using  the 
SYD  at  each  time  gate  for  three  different  object  locations. 
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made  in  the  development  of  sensors  and  inversion  and  classification 
algorithms.  The  first  test  (Shubitidze  et  al.,  2010c;  Fernandez  et  al., 
2010),  held  in  2007  at  the  former  Camp  Sibert  in  Alabama,  was 
relatively  straightforward:  It  was  required  to  discriminate  between 
large,  intact  4. 2"  mortar  rounds  and  smaller  explosion  byproducts 
and  assorted  debris.  The  data  were  taken  by  two  first-generation 
sensors,  the  EM-61  and  EM-63,  both  developed  by  Geonics, 
Inc.  There  were  150  cells  plus  66  for  calibration,  each  one  well 
separated  from  the  others  and  containing  only  one  target. 

The  second  test  (Shamatava  et  ah,  2010;  Shubitidze  et  ah,  2010a) 
was  held  in  2009  near  San  Luis  Obispo  in  California.  This  new  test 
better  reflected  the  wide  range  of  dangerous  targets  usually  present 
in  the  field:  there  were  samples  of  60-mm  mortars,  2.36"  rockets 
(whole  and  fragments),  81 -mm  projectiles,  and  4.2"  mortars;  three 
additional  munition  types  were  discovered  during  the  course  of  the 
demonstration.  Each  cell  contained  an  unspecified  number  of  dan¬ 
gerous  and/or  innocuous  targets.  The  participants  were  expected  not 
only  to  identify  the  UXO,  but  also  to  classify  them  by  caliber.  The 
test,  moreover,  featured  a  more  demanding  topography.  Magnet¬ 
ometers  and  commercially  available  first-generation  sensors  were 
used  to  detect  and  flag  anomalies  that  were  then  interrogated  more 
closely  by  state-of-the  art  EMI  instruments  (TEMTADS,  the  Ber¬ 
keley  UXO  Discriminator,  and  the  Geometries  MetalMapper). 
The  test  was  at  a  much  larger  scale,  comprising  over  1000  data  sets 
for  each  instrument. 

The  third  test  (SERDP,  2010),  on  which  we  concentrate  here, 
took  place  at  the  former  Camp  Butner  in  North  Carolina.  Camp  But- 
ner  functioned  as  a  training  ground  for  infantry  divisions  and  artil¬ 
lery  and  engineering  units  during  World  War  II.  A  large  variety  of 
munitions  were  reportedly  used  there;  at  the  test  site,  the  targets  of 
interest  consisted  mostly  of  37-mm  projectiles  (with  and  without  a 
copper  driving  band)  and  105 -mm  howitzer  shells  and  HEAT 
rounds.  The  clutter  items  were  typical  explosion  byproducts,  such 
as  partial  shells  and  fuzes,  along  with  smaller  shrapnel  and  nonord¬ 
nance  metallic  debris.  Many  of  the  fragments  were  similar  in  size 
and  wall  thickness  to  the  37-mm  projectiles  (Andrews  et  al.,  2011). 
The  test  started  with  an  exploratory  cart-based  survey  using  the 
Geonics  EM-61  sensor.  The  cells  in  which  the  EM-61  data  indicated 
the  presence  of  an  anomaly  (with  one  or  more  metallic  targets)  were 
then  subjected  to  cued  examination  with  TEMTADS.  A  total  of 
2291  cells  were  interrogated,  of  which  171  corresponded  to  poten¬ 
tially  hazardous  targets.  The  discrimination  results  obtained  by  our 
group  are  published  elsewhere  (Cazares  et  al.,  2011).  Here,  we  de¬ 
monstrate  only  our  usage  of  JD  to  estimate  numbers  of  targets  and 
illustrate  shape  differences. 

Figure  8  shows  typical  examples  of  the  two  kinds  of  105-mm 
projectiles  found  at  the  site,  along  with  the  time-dependent  MSR 
eigenvalues  extracted  from  the  cells  that  contained  them.  There 
are  three  above-noise  eigenvalues  in  each  case,  revealing  the  pre¬ 
sence  of  one  target.  The  howitzer  shell  of  Figure  8a  (the  large  TOI  in 
the  test- stand  experiment  described  above)  is  buried  55  cm  below 
ground,  whereas  the  HEAT  round  in  Figure  8b  is  62  cm  deep.  The 
signal  from  the  latter  is  weaker  because  of  the  greater  depth,  despite 
the  lower  noise  level  (which,  in  the  case  of  Camp  Butner,  varied 
substantially  from  cell  to  cell).  The  slow  decay  of  the  eigenvalues 
indicates  that  these  are  large  targets  with  substantial  metal  content 
(Shubitidze  et  al.,  2010b).  Two  of  the  eigenvalues  are  again  seen  to 
be  roughly  parallel,  and  in  the  case  of  the  shell,  they  overlap  each 
other.  The  behavior  of  the  other  (nonparallel)  above-threshold 


eigenvalue  allows  preliminary  classification  of  the  target.  Figure  8c 
shows  the  eigenvalue  decay  curves  (computed  using  both  JD  and 
constrained  nonlinear  optimization,  as  discussed  in  Appendix  B) 
for  the  HEAT  round  as  measured  by  TEMTADS  in  the  APG  test 
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Figure  7.  MSR  eigenvalues  of  the  test-stand  measurements  as  the 
munitions  of  Figure  3  are  placed  together  (a)  and  then  joined  by  the 
sphere  (b)  and  by  the  sphere  and  spheroid  (c)  of  Figure  4. 
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Figure  8.  (a)  MSR  eigenvalues  for  Camp  Butner  Cell  42,  which 
contained  a  105 -mm  HE  shell  buried  at  a  depth  of  55  cm. 
(b)  MSR  eigenvalues  for  Camp  Butner  Cell  62,  containing  a 
105 -mm  HEAT  round  buried  at  62  cm.  (c)  MSR  eigenvalues  for 
a  105-mm  HEAT  round  placed  53  cm  below  TEMTADS  in  the  test 
stand,  computed  using  JD  (solid  lines)  and  nonlinear  constrained 
optimization  (dots). 


stand  of  the  previous  section;  the  munition  was  53  cm  below  the 
sensor  center  in  that  case.  The  time  decay  in  Figure  8b  resembles 
this  measurement  more  than  it  does  the  other  cell. 

Typical  instances  of  the  other  Camp  Butner  TOI,  the  two  different 
kinds  of  37-mm  projectiles,  are  shown  along  with  their  time-depen¬ 
dent  eigenvalues  in  Figure  9.  The  signal  in  Figure  9a  is  strong  and 
distinct,  reflecting  the  shallowness  of  the  target  (it  is  only  8  cm 
deep).  The  copper  driving  band,  with  its  high  conductivity,  influ¬ 
ences  the  signal  by  retarding  its  decay.  The  signal  in  Figure  9b 
is  much  weaker,  as  the  target  is  at  a  larger  depth,  but  it  has  three 
time-dependent  eigenvalues  that  are  distinct  from  the  noise,  whose 
level  is  similar  to  that  of  Figure  9a,  through  the  whole  decay. 

Figure  10a  depicts  the  result  of  running  JD  on  a  cell  that  con¬ 
tained  a  fuze  and  a  small  detonation  fragment,  both  located 
16  cm  below  the  surface,  whereas  Figure  10b  does  the  same  for 
a  single  fuze  buried  almost  twice  as  deep  at  28  cm.  In  both  cases, 
only  three  eigenvalue  decay  curves  are  above  the  noise  threshold, 
indicating  that  the  targets  in  each  cell  can  be  represented  by  a  single 
source.  The  JD  analysis  in  Figure  10a  case  fails  to  identify  the  small 


Time  (ms) 


Time  (ms) 

Figure  9.  (a)  MSR  eigenvalues  for  Camp  Butner  Cell  124,  with  one 
kind  of  37-mm  projectile  at  8  cm  below  ground,  (b)  MSR  eigen¬ 
values  for  Camp  Butner  Cell  192,  which  contained  the  other  kind 
of  37-mm  projectile,  this  one  at  a  depth  of  20  cm. 
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fragment,  which  was  located  close  to  the  fuze.  In  Figure  10b,  there 
are  still  three  dominant  eigenvalues,  but  the  low  signal-to-noise 
ratio  makes  them  barely  discernible.  The  shape  information  given 
by  JD  was  crucial  in  helping  to  identify  this  target  correctly.  Here,  as 
in  all  real-world  situations,  the  noise  is  due  simultaneously  to  the 
ground,  which  can  be  magnetically  responsive  or  wet  (or  both),  to 
the  sensor,  and  to  objects  too  small  or  distant  to  produce  detectable 
individual  signals  but  whose  collective  response  may  be  significant. 
Again,  it  is  seen  that  the  MSR  eigenvalues  mix  positional  informa¬ 
tion  with  intrinsic  response:  small,  shallow  targets  and  large,  deep 
objects  produce  time-dependent  eigenvalues  of  comparable  mag¬ 
nitudes. 

In  Figure  1 1,  we  have  a  large,  shallow  sample  of  “cultural”  (i.e., 
not-ordnance-related)  debris:  a  plowshare  with  depth  of  only  7  cm. 
The  eigenvalues  are  quite  large,  as  expected,  but  decay  very  fast  in 
time,  indicating  that  this  is  not  a  TOI.  The  final  two  figures  corre¬ 
spond  to  multitarget  scenarios.  In  Figure  12a,  we  have  the  extracted 
eigenvalues  corresponding  to  a  set  of  fragments  of  a  155 -mm  muni¬ 
tion  buried  at  a  depth  of  50  cm.  Two  sources  suffice  to  describe  the 


Time  (ms) 

Figure  11.  MSR  eigenvalues  for  Camp  Butner  Cell  14,  with  a  siz¬ 
able  plowshare  buried  at  7  cm. 


Time  (ms)  Time  (ms) 

Figure  10.  (a)  MSR  eigenvalues  for  Camp  Butner  Cell  58,  which  Figure  12.  (a)  MSR  eigenvalues  for  Camp  Butner  Cell  13,  contain- 

contained  two  explosion  byproducts,  a  fuze  and  a  shrapnel  frag-  ing  fragments  of  a  155-mm  projectile,  buried  50  cm  below  the 

ment,  at  a  depth  of  16  cm.  (b)  MSR  eigenvalues  for  Camp  Butner  ground,  (b)  MSR  eigenvalues  for  Camp  Butner  Cell  45,  containing 

Cell  1728,  which  contained  another  fuze  (buried  28  cm  below  numerous  fragments,  buried  at  several  depths,  of  a  munition  that 

ground).  exploded. 
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objects.  Figure  12b  depicts  the  contents  and  data  of  a  veritable  mul¬ 
titarget  cell,  with  several  explosion  byproducts  at  a  wide  assortment 
of  depths.  The  resulting  eigenvalue  map  can  immediately  alert  op¬ 
erators  that  there  may  be  several  potentially  dangerous  targets  in¬ 
terred  below  and  warn  analysts  from  the  outset  that  any  inversion  to 
be  performed  must  involve  a  forward  model  that  assumes  many  tar¬ 
gets.  (Methods  to  transmit  this  information  could  be  incorporated 
into  the  sensor,  and  those  cells  will  always  be  dug  out.) 

The  results  of  the  Camp  Butner  blind  test  were  scored  indepen¬ 
dently  by  the  Institute  for  Defense  Analyses  (Cazares  et  al.,  2011; 
SERDP,  2011).  We  correctly  characterized  all  of  the  2291  anoma¬ 
lies,  each  containing  at  least  one  target,  for  which  data  were 
provided.  All  of  the  171  potentially  dangerous  targets  were  identi¬ 
fied  as  such  and  also  classified  by  caliber,  and  our  algorithms  more¬ 
over  discriminated  between  the  two  kinds  of  37-mm  projectiles  and 
distinguished  the  105 -mm  HEAT  rounds  from  the  other  105 -mm 
shells.  We  asked  for  116  targets  to  be  dug  before  we  were  sure  that 
we  had  identified  all  the  TOI;  in  contrast,  the  best  analysis  using  the 
dipole  model  required  63  extra  digs  to  reach  that  level,  whereas 
other  procedures  needed  many  more  false  alarms  to  be  unearthed 
and  in  some  cases  failed  to  identify  UXO  (Cazares  et  al.,  2011). 
One  of  the  stated  objectives  of  the  Camp  Butner  test  was  to  “pave 
[ . . .  ]  the  way  for  reduced  costs  and  an  accelerated  timeline  to  re¬ 
mediate  munitions-contaminated  sites  throughout  the  nation” 
(SERDP,  2010).  The  success  of  our  procedure  gives  an  idea  of 
the  progress  that  has  been  made  toward  fulfilling  this  aim. 


CONCLUSION 

In  this  paper,  we  applied  a  procedure  based  on  joint  diagonaliza- 
tion  that  provides  a  fast  and  reliable  inversion-free  estimate  of  the 
number  of  targets  buried  under  an  EMI  sensor.  The  number  of  tar¬ 
gets  can  then  be  input  into  multitarget  inversion  procedures,  which 
work  much  faster  and  more  reliably  when  given  this  information. 
The  eigenvalue  decay  curves  found  by  JD  have  amplitudes  that  de¬ 
pend  on  the  location  and  orientation  of  the  targets,  but  their  shapes, 
including  the  appearance  of  parallel  curves  for  cylindrically  sym¬ 
metric  targets,  are  independent  of  those  particulars.  We  reiterate  that 
the  JD  procedure  cannot  replace  inversion,  which  is  still  necessary 
because  the  depths  at  which  the  targets  are  buried  must  be  known  so 
they  can  be  dug  out  safely;  better  knowledge  of  the  locations,  more¬ 
over,  allows  better  determination  of  the  intrinsic  signatures  of  the 
targets  and  consequently  more  robust  and  reliable  identification  and 
classification. 

The  JD  technique  has  several  other  potential  uses,  which  we  have 
started  to  pursue.  It  could  be  used,  for  example,  to  remove  noise 
from  data  by  removing  the  smaller  extracted  eigenvalues  and  their 
corresponding  eigenvectors.  At  the  other  end,  removing  the  larger 
eigenvalues  and  eigenvectors  could  make  it  possible  to  resolve 
smaller  or  deeply  buried  objects.  For  these  applications,  it  is  critical 
to  associate  each  eigenvalue  with  its  eigenvector  unambiguously.  In 
addition,  a  JD  analysis  performed  in  the  field  can  immediately  alert 
data  collectors  about  the  noise  level,  which  can  then  be  adjusted  by 
having  TEMTADS  sample  signals  for  a  longer  time. 

Our  aim  has  been  to  introduce  the  JD  algorithm  and  its  applica¬ 
tion  to  a  particular  sensor,  the  TEMTADS  array,  whose  design 
features  make  the  JD  implementation  particularly  transparent. 
We  have  adapted  the  JD  technique  to  other  EMI  instruments  — 
the  Geometries  MetalMapper  and  a  2  x  2  portable  version  of 


TEMTADS  —  and,  in  the  future,  we  will  present  the  results 
obtained  for  blind  tests  carried  out  with  these  tools. 
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APPENDIX  A 

SIMPLIFIED  MODEL  OF  THE  MSR  MATRIX 


We  consider  an  idealized  version  of  TEMTADS  consisting  of  a 
5x5  array  of  vertical  dipoles  of  unit  moment  z  —  corresponding 
to  a  combination  of  unit  current  and  unit  area  —  each  of  which  at  a 
given  time  gate  transmits  a  primary  field 


Hpr  =  — [3(z  •  R)R  -  zR2] 

4kR5  K  ’  J 

_  3ZX  „  3 ZY  „  3Z2  -  R2  „ 

~  4kR5  X  +  4nR5  y  +  4 kR5  Z’ 


(A-l) 


where  R  =  Xx  +  Yy  +  Zz  =  r  —  rt,  r  is  an  observation  point,  and  rt 
is  the  location  of  the  transmitter.  Each  TEMTADS  dipole  receives 
from  a  buried  source  a  response  Hsc  whose  z-component  reads 
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[3(z  •  R)R  —  z R2]  ■  m, 
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where  now  R  =  ly  -  r  with  rr  the  location  of  the  receiver.  Note  that 
the  expressions  in  square  brackets  in  equations  A-l  and  A-2  coin¬ 
cide  if  and  only  if  the  (infinitesimal)  transmitters  and  receivers  are 
colocated,  as  is  the  case  with  TEMTADS.  The  moment  m  induced 
in  the  responding  source  is  given  as  in  equation  4.  The  MSR  matrix 
for  the  transmitter/receiver  array  is  constructed  by  assembling  Nt 
samples  of  equation  A-l  for  the  transmitted  field  and  Nr  samples 
of  the  received  response  A-2,  thus  forming  the  Green  tensors 
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The  measured  secondary  signal  given  in  equation  1  is 

S  =  Gsc(UAUr)(Gpr)r  =  (GscU)A(GprU)r,  (A-5) 


where  U  and  A  are  3  X  3  matrices. 

We  can  use  this  simplified  model  to  provide  a  plausibility  argu¬ 
ment  for  one  of  our  findings:  Targets  with  cylindrical  or  spherical 
symmetry  have  a  pair  of  eigenvalues  whose  log-log  decay  curves 
look  parallel  (i.e.,  there  is  a  constant  ratio  between  the  two).  Take 
Gpr  and  Gsc  to  be  identical,  as  the  geometry  above  allows,  and 
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consider  a  2  x  2  sensor  array  such  that  the  Tx/Rx  pairs  are  located  at 
17  =  <f(±x  ±  y),  where  €  =  0.20  cm  in  the  case  of  TEMTADS. 
Suppose  the  target  is  a  dipole  buried  at  depth  z0  below  the  center 
of  the  array  and  characterized  by  the  cylindrically  symmetric  polar¬ 
izability  tensor  A  =  /?diag(l,  1,2).  Ordering  the  Tx/Rx  dipoles 
clockwise  from  the  top  left  of  the  array,  we  obtain  the  MSR  matrix 

-11  1-1 

11-1-1, 

K  K  K  K 

(A-6) 


where  I  is  the  identity  matrix.  We  accomplish  this  by  making  re¬ 
peated  Givens- Jacobi  similarity  transformations  designed  to  gradu¬ 
ally  accumulate  the  “content”  of  the  matrices  on  their  diagonals 
until  a  certain  tolerance  level  is  reached.  The  transformations  are 
of  the  form  A {tq)  ->  A '{tq)  =  \ rsA(tq)Yjs,  with  the  matrix  Yrs 
being  the  identity  but  with  the  four  elements  Vrr,  Vrs,  Vsr,  and 
Vss  replaced  by  the  2D  rotation  array 

cos  cp  sin^J  tan2(prs  =  - 

sin  (j)rs  cos  (f)rs  J  +  yj f2rs  +  n2rs 

(B-2) 
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where  we  have  introduced  the  dimensionless  parameter 


_  3zp  -  fl2  =  2  /zo  _ 
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which  implies 


(A-7) 


If  the  target  is  a  horizontal  dipole,  the  polarizability  tensor  is 
UAUr  =  /?diag(2, 1,1)  and  the  MSR  matrix  A-6  becomes 
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(1  +  2)  +  /e2  _ 
(A- 8) 


The  eigenvalues  of  this  4x4  matrix  are  0,  4,  4/c2,  and  42.  One  ei¬ 
genvalue  vanishes,  confirming  that  there  is  only  one  target.  The 
nonzero  eigenvalues  are  in  general  not  degenerate,  even  for  highly 
symmetric  configurations.  (The  eigenvalues  coincide  for  the  special 
case  z0  =  2if.)  On  the  other  hand,  there  is  a  constant  ratio  k 2  be¬ 
tween  two  of  the  eigenvalues,  allowing  extraction  of  the  extrinsic 
depth  zo  through  the  definition  in  equation  A-7.  (Cylindrical  sym¬ 
metry  is  important  here;  if  A  =  /?diag(l  ,^,2),  the  MSR  eigenva¬ 
lues  become  0,  4/c2,  4^,  and  42,  with  all  significant  ratios  now 
time-dependent.)  This  undertaking  quickly  becomes  nontrivial: 
for  example,  when  the  target  is  vertical,  the  MSR  eigenvalues 
are  0,  4,  4,  and  42/c2;  the  intrinsic  and  extrinsic  features  are  coupled. 
The  analysis  is  further  complicated  if  the  target  is  displaced  from  the 
center  of  the  Tx/Rx  array  or  has  another  orientation,  or  if  the  Tx/Rx 
array  is  larger  than  2x2. 

APPENDIX  B 

ALGORITHM  FOR  JOINT  DIAGONALIZATION 

The  joint  diagonalization  algorithm  we  use  (Flury  and  Gautschi, 
1986;  Cardoso  and  Souloumiac,  1996;  Belouchrani  et  al.,  1997)  is  a 
generalization  of  Jacobi’s  procedure  to  find  the  eigenvalues  of  a 
single  matrix.  Formally,  we  solve  the  optimization  problem 

nnn  5  E  ([VA(f9)Vr](/)2 ,  ^ 

s.t.  VrV  =  I 


Ylys  ^  ^jWP'rritq)  ^ss{tq)\  \^rs{^q) 

q 

(B-3) 

f rs  2^^[tZrr(/g)  Clss{tq)\\cLrs{^q)  (B-4) 
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The  indices  are  swept  systematically  and  the  procedure  is  re¬ 
peated  until  convergence  is  reached.  The  computational  burden 
is  equivalent  to  that  of  diagonalizing  the  matrices  one  by  one. 
The  resulting  eigenvalues  and  eigenvectors  are  all  real  because 
all  the  MSR  matrices  are  symmetric. 

We  have  also  solved  the  nonlinear  constrained  optimization  pro¬ 
blem  B-l  directly  using  sequential  quadratic  programming  as  im¬ 
plemented  in  the  function  fmincon  of  the  Matlab  Optimization 
Toolbox  (The  MathWorks,  2006).  This  procedure  unsurprisingly 
takes  a  longer  time  to  run  than  the  Givens-Jacobi  algorithm 
sketched  above.  The  two  procedures  gave  equivalent  results  in 
all  cases  we  tried;  Figure  8c  shows  an  example. 
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